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We use quantum Monte Carlo (stochastic series expansion) and finite-size scaling to study the 
quantum critical points of two 5=1/2 Heisenberg antiferromagnets in two dimensions: a bilayer 
and a Kondo-lattice-like system (incomplete bilayer), each with intra- and inter-plane couplings J 
and J± . We discuss the ground-state finite-size scaling properties of three different quantities — the 
Binder moment ratio, the spin stiffness, and the long-wavelength magnetic susceptibility — which we 
use to extract the critical value of the coupling ratio g — J±/ J. The individual estimates of g c 
are consistent provided that subleading finite-size corrections are properly taken into account. In 
the case of the complete bilayer, the Binder ratio leads to the most precise estimate of the critical 
coupling, although the subleading finite-size corrections to the stiffness are considerably smaller. 
For the incomplete bilayer, the subleading corrections to the stiffness are extremely small, and this 
quantity then gives the best estimate of the critical point. Contrary to predictions, we do not find a 
universal prefactor of the ~ 1/L spin stiffness scaling at the critical point, whereas the Binder ratio 
is consistent with a universal value. Our results for the critical coupling ratios are g c = 2.52181(3) 
(full bilayer) and g c — 1.38882(2) (incomplete bilayer), which represent improvements of two orders 
of magnitude relative to the previous best estimates. For the correlation length exponent we obtain 
v = 0.7106(9), consistent with the expected 3D Heisenberg universality class. 

PACS numbers: 75.10.Jm, 75.10.-b, 75.40.Cx, 75.40.Mg 



I. INTRODUCTION 

The two-dimensional (2D) S —1/2 Heisenberg antifer- 
romagnet has received considerable attention in the past 
two decades because of its close relation to the Cu02 lay- 
ers of the cuprate superconductors^ Other, even better 
realizations of this model system have been discovered as 
welli^ The properties of the single-layer Heisenberg model 
have been thoroughly studied using both analytical and 
numerical methods, and there is very good agreement 
with experiments, e.g., for the temperature dependence 
of the spin correlation length&i in La2Cu04 (measured 
using neutron scattering) and NMR relaxation rates£ 

Mapping the lattice Heisenberg model onto a contin- 
uum field theory yields the (2+l)-dimensional nonlinear 
sigma-modelj^i the coupling constant g of which con- 
trols the transition from Neel order to quantum disorder 
at temperature T = (quantum phase transition 7 ). This 
transition is in the universality class of the finite-T tran- 
sition of the 3D classical Heisenberg model. 3,8 Having an 
ordered ground stated the 2D square-lattice S = 1/2 
Heisenberg model corresponds to g < g c . Even so, there 
is some influence from the critical point, because a quan- 
tum phase transition is also associated with universal 
quantum critical scaling at finite temperature, in an ex- 
tended (g, T) regime where temperature is the dominant 
energy scaled The energy scales characterizing the or- 
dered and disordered phases — the spin stiffness and the 
singlet-triplet gap, respectively — vanish continuously as 
g — > g c , and hence the quantum critical regime fans out 
from the point (g = g c , T = 0). 

A quantum phase transition of the type described by 
the nonlinear sigma-model can be realized in the Heisen- 



berg antiferromagnet by introducing a pattern of two (or 
more) different coupling strengths in a way that favors 
singlet formation on dimers (or larger units of an even 
number of spins) AAfi This leads to an order-disorder tran- 
sition at some critical coupling ratio. Models of this kind, 
e.g., a bilayer where dimers form across the layers ^LiSii^ 
single layers with various dimer patterns^ or a regu- 
larly depleted system where singlets form on rings of 
four or eight spins^ have been extensively studied using 
quantum Monte Carlo simulations in order to confirm 
the expected universality class and to test very detailed 
predictions** of the finite-T quantum critical behavior of 
various quantities. The predicted universal behavior was 
confirmed at low temperatureii 3 ^^^^ The simulations 
also served to establish the onset of nonuniversal lattice 
effects at higher temperature and the nature of the cross- 
over— to the low-temperature renormalized classical or 
quantum-disordered regimes away from the critical point. 

In this paper we study the critical points of two differ- 
ent S = 1/2 Heisenberg models: a symmetric bilayer and 
a Kondo-lattice-like system in which there are no intra- 
plane couplings in one of the layers. We will refer to these 
systems as the full and incomplete bilayers, respectively; 
see Fig. ^ Their Hamiltonians {H\ for the full bilayer 
and H2 for the incomplete bilayer) are given by 

Hi — J (Su ■ Sij + S 2 i • S%j) + J± Sh • S 2 i, (1) 

(i,3) 1 

H 2 = J ^ Sh • Sij + Jj_ ^ Su ■ S 2 i- (2) 

(i,3) 

Here, S a ,i is a spin-1/2 operator at site i of layer a 
(a = 1,2), and denotes a pair of nearest-neighbor 
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(a) 

FIG. 1: The arrangement of spin interactions in the (a) full 
and (b) incomplete bilayers. There are two different cou- 
plings: J (intra-plane) and J± (inter-plane), as indicated. 



sites on the square lattice of L x L sites with periodic 
boundary conditions. Both coupling constants are anti- 
ferromagnetic (J, J± > 0). In order to avoid any frustra- 
tion we consider only even L. As the ratio g = Jj_/J is 
increased, there is a tendency to form inter-plane near- 
neighbor singlets, which at some g — g c leads to the 
opening of a spin gap and destruction of the long-range 
Neel order present for g < g c (in the limit g — > oo the 
ground state is a singlet product, which clearly cannot 
support long-range spin-spin correlations). This is the 
transition we investigate in detail in this paper. 

Our purpose in studying these models is two-fold. 
First, we would like to determinate the locations of 
the critical points to much higher accuracy than they 
are currently known. The best estimates to date are 
g c = 2.525(2) (full bilayer)A& and g c = 1.393(8) (incom- 
plete bilayer )Si (see also Ref . Il7j) . The statistical accu- 
racies here are quite modest compared to results for the 
standard classical critical points (e.g., the 3D Heisenberg 
modeliSii*-) . It would be useful to increase the preci- 
sion so that studies of various aspects of finite- T quan- 
tum criticality (e.g., interesting properties of isolated im- 
purities in a critical host systemSLSSiSiSi^) could be 
studied numerically at low T very close to the criti- 
cal point (minimizing the effects of the eventual cross- 
over to the renormalized-classical or quantum-diosrdered 
regime). Second, we wish to compare several different 
ways of extracting the critical coupling, in order to gain 
additional confidence in the results and to provide guid- 



ance for studies of other quantum critical points. The 
reason for choosing the particular bilayer lattices of Fig.Q] 
over other 2D Heisenberg systems undergoing the same 
type of transition is that they do not break any in-plane 
symmetries of the square lattice. 

We have carried out finite-size scaling of low- 
temperature (T — > converged) QMC results for three 
different quantities: the spin stiffness, the Binder cumu- 
lant ratio, and the long-wavelength (uniform) magnetic 
susceptibility. We use our recently proposed method for 
including subleading finite-size corrections^ Although 
this necessitates nonlinear fits with a relatively large 
number of independent parameters, we believe that this 
is necessary in order to obtain completely unbiased re- 
sults. Our final results for the critical couplings are 
g c = 2.52181(3) for the full bilayer and g c = 1.38882(2) 
for the incomplete bilayer, i.e., the statistical precision 
is improved by approximately two orders of magnitude 
relative to the previous estimates. Our fitting proce- 
dure also gives the correlation-length exponent v, but 
because of the multi-parameter fits and the relatively 
modest lattices sizes (L up to 42), its statistical precision 
is not quite as high (the error bars are roughly twice as 
large) as that of recent classical Monte Carlo simulations 
of the 3D Heisenberg model ML Nevertheless, our result, 
v = 0.7106(9), is fully consistent with the presently most 
accurate value of this exponent, v — 0.7112(5)^£ 

We also discuss the universality of the Binder moment 
ratio and the spin stiffness at the critical point. In the 
former case, we point out a difference relative to classical 
systems in how the order-parameter moments are defined 
and calculated, and in the latter case we do not find the 
universality that has been suggested for the prefactor of 
the ~ 1/L scaling^ (i.e., we obtain different prefactors 
for the two different bilayer models). 

The rest of the paper is organized as follows. In Scc.lTTl 
we discuss the quantities that we have calculated and 
their QMC (Stochastic Series Expansion, SSE) estima- 
tors, as well as their expected finite-size scaling forms. 
In Sec. IIIII we first briefly review our approach to deal 
with subleading finite-size corrections and then present 
the results of the analysis. We give a brief summary and 
conclusions in Sec. IIVI 



II. CALCULATED OBSERVABLES AND THEIR 
CRITICAL SCALING PROPERTIES 

We have used the SSE QMC method with operator- 
loop updates^ This approach is based on sampling di- 
agonal matrix elements of the power series expansion of 
exp(—[3H), where (3 is the inverse temperature. We use 
L x L x 2 lattices with periodic boundary conditions in 
the x and y directions, with even L up to 42. In order 
to ensure convergence of all calculated quantities to their 
ground-state values, we carried out simulations at inverse 
temperatures (3 n = 2™ with integer n taken large enough 
so that results for f3 n and /3 n -i agree within statistical 
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FIG. 2: (Color online) Convergence as a function of inverse 
temperature /3 — 2™ of the squared sublattice magnetization 
to its ground state value for three different lattice sizes. 



errors. Examples of the convergence are shown in Fig. [21 



A. Binder Moment Ratios Qk 

The magnetic moment ratios Qk introduced by 
Binder— have the very useful property of being universal 
at the critical point, because all nonuniversal scale fac- 
tors cancel out along with the length dependence. This 
follows from the finite-size scaling hypothesis for the or- 
dered moment (here the staggered magnetization). The 
k th power of the staggered magnetization scales as 



(|m| fe ) L = L-W v M k {tlM v ), 



(3) 



where L is the linear system length, (3 and v are critical 
indices in their standard notation, and t is the reduced 
coupling constant, which we define here in terms of the 
coupling ratio g as t — (g — g c )/9c- Mk{x) are the scaling 
functions. Consequently, the moment ratios 



,2fc\ 



Q k (t,L) = 



(4) 



are dimensionless scaling functions. At the critical point, 
Qk(0, oo ) are universal constants. 

We have computed the first two Binder ratios, which 
we define as 



(m 2 ) _ 3((m z ) 2 ) 
~ 2(|m*|) 2 ' 
5((m z ) 4 ) 



(Mr 

(m 4 ) 



(m 2 r 3((m*) 2 ) z 



(5) 
(6) 



where m z is the z-component of the staggered magneti- 
zation operator, 



= ^ = mcos(e). (7) 



Of 




FIG. 3: Binder ratio Qi for different system sizes vs the 
coupling ratio for the full (a) and incomplete (b) bilayers. 
Results for even L from 8 to 42 are shown (all except for 
L — 22, 26, 34, 38). The slope of the curves increases with L. 



Here, N = 2L 2 is the number of lattice sites. Since the 
0(3) spin-rotational symmetry is not broken in the sim- 
ulations (i.e., an average over all angles is obtained) 
we have included the appropriate factors to compensate 
for the rotational averaging of m z in Eqs. (JSJ and ©. 

In Fig. we show the ratio Q2 for both bilayer systems 
as a function of g for lattices of different linear length 
L. One can clearly see the curve crossings, indicating 
a quantum critical point, but it is apparent that there 
are sizable corrections to their location. We will analyze 
these crossing-point shifts in Sec. III. 

Note that Q2 at the crossing point is approximately the 
same for both models, in accord with an expected uni- 
versal value as L — > 00. Simulations done on 3D classical 
Heisenberg modelsi^^iSi 2 ^ gave a universal value in the 
range 1.35-1.40, i.e., substantially larger than what we 
see in Fig. (clearly these values have not yet completely 
converged to their infinite-size Qi , and the trend is for the 
crossing value to increase with L, but we will show that 
they converge to Qi ~ 1.29). This disagreement with 
the classical value is easily accounted for by considering 
the way the sublattice magnetization is defined and com- 
puted in a quantum system: Although the 2D system 
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formally is mapped onto a 3D classical model, (|m z | fc ) 
is an equal-time expectation value, which in the simu- 
lations is averaged over the third (imaginary time) di- 
rection. This corresponds directly to taking expectation 
values over individual layers in a 3D classical model. We 
are not aware of any such calculation and hence cannot 
compare directly with the corresponding classical univer- 
sal value. Nevertheless, as we will show in greater detail 
in Sec. Ill, the crossing Q 2 values for both our systems 
are fully consistent with each other and hence support 
universality. 

B. Spin Stiffness p s 

In continuum field theory language, a stiffness p is de- 
fined in terms of the increase in free-energy density / as 
a boundary-condition twist $ is imposed on the order- 
parameter field 9: 

8f s {t,L) = \p{Vef = \p{$/Lf. (8) 
The prefactor is the stiffness constant, 



At a quantum critical point, it should scale as^S 

p~L 2 - d -\ (10) 

where d is the dimensionality and z is the dynamic critical 
exponent. Wallin et al. have argued that pL 2 ~ d ~ z is a 
universal number in two dimensions^ 

For the Heisenberg model, the spin stiffness p s is deter- 
mined by imposing a twist directly in the Hamiltonian, 
modifying the spin-spin interactions in one of the lattice 
directions: Sj • Sj — * Si ■ R($/Z)Sj, where R rotates the 
spin operator about an appropriately chosen axis^ In 
SSE simulations, like in path integrals^ the stiffness is 
directly obtained without explicitly imposing a twist, as 
the second derivative of the energy with respect to the 
twist at $ = 0. This leads to an estimator in terms of 
winding number fluctuations^ 

p s = \(w 2 x + w 2 y )/P, (11) 

where the winding numbers are 

w a = (N+ -N~)/L, (a = x,y). (12) 

Here 7V + (N~) is the number of operators S~ (S~ Sj~) 
in the sampled terms of the power series expansion, with 
i,j two nearest-neighbor sites oriented along the lattice 
a (x or y) axis. The definition 1|11[> corresponds to the 
stiffness per unit cell of the bilayer models. 

In the case of the bilayer models we have d = 2 and 
expect z — 1, and hence the scaling i|l(J|) becomes p s ~ 
L _1 . The quantity p s L should thus be size-independent 
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FIG. 4: The spin stiffness multiplied by the system length 
L at the critical point. The system sizes are the same as in 
Fig. |21 The slope of the curves increases with L. 

at the critical point, and also in this case one can expect 
curves for different L to cross each other. Such crossings 
have previously been used to approximately locate the 
critical point of the full bilayer. 32 

Fig. 01 shows our SSE results for p s L versus g for dif- 
ferent lattice sizes. Again, one can see that the crossing 
points move as L is increased, but, interestingly, much 
less so for the incomplete than the complete bilayer. 
These results do not immediately support a universal Lp s 
at g c ; a careful finite-size scaling analysis does not change 
this conclusion, as we will see in Sec. III. 

C. Uniform Susceptibility \{q ~* 0) 

The temperature dependence of the uniform magnetic 
susceptibility is an often-used indicator of quantum crit- 
icality. Exactly at g — g c , the general asymptotic scaling 
form isiS 

Xu(T) ~ T d ' z -\ (13) 

This has been numerically confirmed at low T in previ- 
ous QMC simulationsiiiii and series expansions 3 ^ of the 
bilayer and other critical Heisenberg systems. 35 Here we 
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will consider the corresponding finite-size scaling behav- 
ior, which we obtain by substituting the finite-T quantum 
critical correlation length, £ ~ y-i/z^g^ m ^ ne tempera- 
ture dependence l|13[) and then substitute L for £, giving 
X„ ~ L z ~ d . However, we apparently have a problem here 
since the uniform susceptibility \u — (3{{M Z ) 2 ) vanishes 
as T — > 0, due to the conserved magnetization M z and 
the singlet ground state. In order to carry out finite-size 
scaling, we therefore consider the long-wavelength limit 
of the wave- vector dependent susceptibility x(q)> which 
we obtain in practice by taking q = 2tt/L. Thus we will 
test the finite-size scaling form 

-d 



X(q - 0) = x(2tt/Z) ~ L z 



(14) 



The static spin-spin susceptibility in real space is given 
by the Kubo integral 

X (k,l)= [ dr(S* k (T)S?(0)), (15) 
Jo 

which in SSE simulations is obtained in terms of 
spins in the states propagated by the sampled operator 
sequences^ 



X(k,l) 




^S z k [p]St[p] 
p=o j 



(16) 



n is the number of hamiltonian operators in the sampled 
sequences and the index p refers to the state obtained 
after p operators have acted. The Fourier transform that 
we are interested in is 

x(q) = ^£ eiq ' (rh ~ ri) x(M), ( 17 ) 

k,l 

and since our models are symmetric with respect to a 90° 
rotation, we take the long-wavelength susceptibility as 



Xu = x(« ^ o) = ~ xi'fA 0) + x (o, ¥, 0) 



(18) 



We again consider the form leading to curve crossings 
at the critical point, i.e., we plot XuL, which should be 
size- independent at g c . Figure [5] shows the data that we 
will anayze more carefully in the next section. Again we 
observe crossing points, which shift significantly with L. 



III. DATA ANALYSIS 

We first discuss here a rough determination of the 
critical coupling ratios of the two models, studying the 
asymptotic behavior of the crossing points discussed 
above. This will also serve as a first confirmation of mu- 
tual consistency of the leading scaling forms for the three 
different quantities under consideration. We then analyze 
the results in greater detail using a finite-size scaling hy- 
pothesis including subleading corrections. 
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FIG. 5: The long- wavelength susceptibility versus the cou- 
pling ratio for (a) the full bilayer and (b) the incomplete bi- 
layer. The system sizes are those listed in Fig. The slope 
of the curves again increases with L. 



A. Critical Coupling from Crossing Points 

We use the data presented in the previous section to 
extract the intersection points of fixed- L curves for sys- 
tem sizes L and 2L (other size ratios give similar results) . 
Our simulations have been performed on a rather dense 
grid of g-points, and we can therefore reliably obtain the 
intersection points using fits of straight lines or second- 
order polynomials to interpolate between the data points. 
In Fig. we plot the results versus the inverse system 
size. For both models, the crossing points drift toward 
a common critical coupling in the L — > oo limit, thus 
confirming the scaling laws discussed in the prevcious 
section. 

For both models, especially the incomplete bilayer, the 
spin stiffness crossing point exhibits the most rapid con- 
vergence (i.e., the weakest subleading corrections). It 
and the susceptibility converge from above, while the 
Binder ratios converge from below. We can hence bracket 
g c using these results. However, a much more pre- 
cise bracketing can be obtained from the spin stiffness 
curves alone, noting that they become very flat as L 
grows. With the curvature decreasing with increasing L, 
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FIG. 6: (Color online) Convergence vs the inverse lattice size 
of intersection points of curves for L and 2L, for the full (a) 
and incomplete (b) bilayers. The error bars are are much 
smaller than the symbols. The circles at 1/L — indicate the 
critical couplings from the careful finite-size scaling analysis 
carried out in Sec. Ill B. 

a straight-line extrapolation using a few large- L points 
(we use four) should give a lower bound for g c , while 
the crossing point for the largest L should be an up- 
per bound. The critical couplings extracted this way are 
g c e (2.5205,2.5232) and g c G (1.38870,1.38895) for the 
full and incomplete bilayers, respectively. These values 
are fully consistent with the best previous estimates, dis- 
cussed in Sec. I, but the precision is significantly higher. 
The more rigorous data analysis discussed below will fur- 
ther improve on the accuracy. 

Naively, one might expect that the asymptotic ap- 
proach of the crossing points to the critical coupling 
should be given by the correlation-length exponent v, 
as Across = ffc + aL~ x / u , as is the case for fixed-size es- 
timates of the critical coupling (or the critical tempera- 
ture), such as the location of the maximum of the order- 
parameter susceptibility or the specific heat (in the case 
of finite- T transitions). However, a crossing point can- 
not be regarded as a conventional fixed-L definition of 
g c , since two system sizes are involved and there can be 
cancelations of a leading behavior defined in terms of the 
individual lattices. Thus, we would in general expect the 
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FIG. 7: (Color online) The blue crosses (x) show the input 
points of the optimization procedure. The red plusses (+) 
show the output. The dashed black line indicates the magni- 
fied region correspondig to the bottom-right plot in Fig. |5] 

crossing points to converge faster than L~ x l v . Binder 
has discussed the corrections to the cumulant crossing- 
points and in a recent article^ we have presented a 
different way of analyzing crossing points in general (i.e., 
not only for Binder ratios but also for other quantities 
that are size- independent at g c ) which takes subleading 
finite-size corrections into account explicitly. There we 
also showed some results for the spin stiffness crossings 
of the full bilayer model. Here we will not analyze the 
crossing points in any greater detail, but instead consider 
the scaling of the full fixed-L curves shown in the figures 
of Sec. II. Such "data collapse" makes better use of the 
full range of simulation results and can also be carried 
out with subleading corrections taken into account^ 

B. Finite-Size Scaling with Subleading Corrections 

The scaling ansatz typically used to analyze finite-size 
data A(t, L) for a quantity A at reduced coupling t = 
(g — g c )/ g c on a lattice of length L is 

A(t,L)=L K /»f A (tL 1 /»), (19) 

where k is a critical exponent which depends on the quan- 
tity A, i.e., A(t,L = oo) ~ t~ K . This form can be used 
to collapse data in a neighborhood of t — 0, by graphing 
A(t, L)L K I V versus tL -1 /" , adjusting g c , k, and v to ob- 
tain the tightest collapse of the data onto a single curve. 

In Ref. |2fj we started from renormalization group the- 
ory and derived an extension to (|19f) that includes both 
"shift" and "renormalization" corrections: 

A(t, L) = L K ' V {1 + cL-")g A (tL l l v + dL^"). (20) 



7 



1.3890 



1.3889 - 



1.3887 




0.700 



1.3890 



1.3887 



t r 

0.705 0.710 



0.715 




0.700 



0.705 



0.710 



0.715 



2.5227 



2.5225 - 



2.5223 - 



2.5221 - 



2.5219 - 



2.5217 -| 1 1 r 

0.700 0.710 0.715 0.720 0.725 



2.52190 



2.52185 - 




g c 2.52180 



2.52175 - 



2.52170 




i 1 r 

0.7100 0.7125 0.7150 0.7175 0.7200 



FIG. 8: (Color online) The number density of bootstrapped (f,g c ) solutions is plotted for Q2 and p s in greyscale. The left 
panel shows the incomplete bilayer and the right panel the complete bilayer. The solid red lines show the contour at 1/3 the 
maximum density (this would correspond to 2/3 of the weight under a gaussian distribution). In the top right panel, note a 
two-peak structure which sometimes appears in the nonlinear fits. The dotted lines are guides to the eye. They show the Q2 
contour on the p s plot and vice versa. 



Here, <f> is the subleading irrelevant RG eigenvalue, which 
causes a shift in the critical coupling, to is an effective 
exponent that accounts for corrections due to the inho- 
mogeneous part of the free energy and nonlinearity of the 
scaling fields. The constants c and d are nonuniversal and 
should be regarded as fitting parameters along with the 
leading and subleading exponents. From Eq. l]2Up. we 
see that we can now achieve data collapse by plotting 
A(t, L)L~ K I V '/(l + cL~ u ) versus x = tL l / u + dL-^i" for 
different sizes L. 

To carry out this type of analysis in practice, we note 
that the scaling function g A (x) is well-behaved and can 



be Taylor expanded close to the critical point: 

g A (x) = A(t,L)L- K /»/(l + cL-») (21) 
= q + qix + q 2 x 2 + q 3 x 3 + q A x A H 

For the Heisenberg bilayers, the critical indices k and v 
are expected to be those of the 3D classical Heisenberg 
universality class. In the case of the ratios we are consid- 
ering, k/v = are known integers which we hence do not 
have to adjust. The current best estimate for the cor- 
relation length exponent is v = 0.7112(5))** but in our 
analysis we keep it as a free parameter, along with g c , 
the subleading exponents u> and 4>, the constants c and 
d, and the parameters of the polynomial in l|21|) . This 
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amounts to a large number of fitting parameters, but it 
should be noted that the polynomial expansion of the 
scaling function is essentially just an interpolation of the 
data and hence is highly constrained; the freedom in the 
coefficients qt do not add significant freedom to the other 
paramers of the fit (we use a quartic polynomial). More- 
over, the number of fitting parameters is dwarfed by the 
number of data points to be fit (hundreds or thousands) . 

Nonlinear curve fitting has well-known problems as- 
sociated with the convergence of the parameters to the 
globally optimal fit. In our work we already know rather 
accurate estimates for v and g c , and at the first stage of 
the fits we used those values as initial conditions. Once 
we obtained rough estimates for c, u), d, <f>, we used the fol- 
lowing procedure: Performing bootstrap sampling of the 
raw data, we carried out a large number (typically around 
1000) of fits with initial conditions for all the parameters 
taken at random from inside a "box" in parameter space. 
This box is determined such that the fits converge well, 
but that the variation in starting points is significantly 
larger than the final spread of the resulting parameters. 
We then use the spread among the bootstrap samples to 
calculate statistical errors. This procedure is illustrated 
in Figs. □ and il 

The scaling formula, Eq. J20Jl, is strictly valid only for a 
small range of couplings in the vicinity of g c [although the 
range of validity should be larger than with the leading- 
order form l|19l) ]. The parameters obtained show some 
dependence on the range of data included. In order to 
eliminate as much as possible potential remaining effects 
of further subleading corrections that are not captured by 
our approach, we chose to use a rather narrow window in 
the scaled coupling x = (g — g^L 1 ^ /g c + dL~' t '/ u , so that 
there is no longer any statistically detectable dependence 
on the size of the window. Our final results are based on 
x e [—0.25,0.25]. There are also other subtle issues in 
the fitting procedure, e.g., for a given range of x, differ- 
ent number of data points are available for the different 
lattice sizes, typically leading to relatively smaller sta- 
tistical weight for the larger sizes than the smaller sizes. 
We therefore made sure to include only system sizes suffi- 
ciently large for the extracted parameters not to change 
appreciably when systematically excluding more of the 
smaller lattices. We kept only L > 8 for the results we 
report here. 

We found both the Binder ratio and the spin stiffness 
to be well behaved in the fitting procedures. The result- 
ing collapsed data for these quantities, i.e., their scaling 
functions gA(x) in Eq. I|22|l . are shown in Figs. l9l and 1 101 
We do not discuss results for Qi here as it behaves sim- 
ilarly to Q2 and is statistically strongly correlated to it. 
The long- wavelength susceptibility also exhibits data col- 
lapse, however with substantially larger fluctuations than 
p s and Q2 (the prefactor d appears to be rather large and 
difficult to determine accurately, and thus it is difficult 
to fix the data window x € [x m - m , a; max ] in a meaningful 
way). We have therefore focused on p s and Q2 for the 
final high-precision statistical analysis. 
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FIG. 9: (Color online) Data collapse of the Binder ratio for 
the full (a) and incomplete (b) bilayer. The values of g c and 
v obtained are listed in Table U] 



The final parameters and their statistical errors were 
determined from the bootstrap samples. The distribu- 
tions are illustrated by density plots in Fig. [S] We list 
the values for <? c , v, and the value of the respective quanti- 
ties at the critical point, q c = Q2(g c ), Ps(gc)L, in Tabled 
The critical couplings are here consistent among all the 
fits, and the correlation length exponent is marginally 
consistent (within 2-3 standard deviations). 

As seen in the table, The highest relative precision of 
g c is obtained using Q2 for the full bilayer and p s for the 
incomplete bilayer. The latter can probably be traced to 
very small subleading corrections, as is evident already in 
Fig. For the full bilayer p s also has smaller subleading 
corrections than Q2, but still we obtain a higher accuracy 
in g c using Q2 ■ This shows that small subleading correc- 
tions are not necessarily advantageous; what matters is 
how well those corrections are described by the finite-size 
scaling forms used. This in turn should be determined 
by the extent of corrections of even higher order. 

For the subleading exponents lu and (f> we obtain values 
around unity, except in the case of the p s scaling of the 
incomplete bilayer, which requires larger values. For the 
full bilayer, w = 1.14(3), (j) = 0.8(2) in the Q 2 scaling 
and lu = 1.0(3), 4> — 1-2(2) in the p s scaling. For the 
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FIG. 10: (Color online) Data collapse of the spin stiffness 
ratio for the full (a) and incomplete (b) bilayer. The values 
of g c and v obtained are listed in Table [I] 



incomplete bilayer, lu — 1.0(4), <j> — 1.3(2) in Q 2 and 
u> = 1.9(2), (f> — 1.8(2) in p s . All these subleading expo- 
nents should be interpreted as effective ones, as they are 
to some extent affected by the higher-order corrections 
that we have neglected. The fact that u) and <fi obtained 
from p s of the incomplete bilayer are close to 2 suggests 
that in this case the leading corrections are small and the 
exponents instead reflect predominantly the corrections 
of the next higher order. 



IV. SUMMARY AND CONCLUSIONS 

We have carried out finite-size scaling analyses of high- 
precision stochastic series expansion QMC data for two 
S = 1/2 Heisenberg bilayer models. Using three differ- 
ent quantities, the Binder order parameter moment ratio, 
the spin stiffness, and the long-wavelength magnetic sus- 
ceptibility, we obtained very accurate estimates for the 
critical couplings. We have stressed the importance of 
including subleading corrections in the finite-size scaling 
analysis. All the quantities considered then give mutu- 
ally consistent results for the critical couplings as well as 
for the correlation-length exponent v. We have assumed 



TABLE I: Bootstrap averages and errorbars for the param- 
eters g c ,v,qo (the values of the respective quantities at g c ). 
The top group of values are for the incomplete bilayer, and 
the bottom group for complete bilayer. The indicated error 
bars correspond to one standard deviation of the probability 
distributions obtained in the bootstrap analysis, as explained 
in the text. 
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go 


Qi 


1.38885(5) 


0.708(2) 


1.293(3) 


p s L 


1.38882(2) 


0.705(2) 


0.434(3) 


XuL 


1.388(1) 


0.7(1) 


0.28(2) 


Q 2 


2.52180(3) 


0.715(2) 


1.2858(3) 


p s L 


2.5221(2) 


0.714(5) 


1.13(3) 




2.521(1) 


0.7(1) 


0.12(2) 



that the dynamic exponent z — 1 , and all our results are 
completely consistent with this expectation^ 

The inclusion of two different subleading corrections in 
the data fits implies larger statistical fluctuations com- 
pared to an analysis neglecting subleading corrections or 
taking them into account less completely than we have 
done here. However, because of these procedures, along 
with careful studies of the dependence on the range of 
system sizes and the data window around the critical 
point, we are confident that the remaining errors are 
purely statistical and accurately captured by the error 
bars quoted in Table HJ We also note that the different 
quantities we have considered correspond to averaging 
functions of very different properties of the configurations 
generated in the simulations — the staggered magnetiza- 
tion in the case of the Binder ratio, the winding number 
in the case of the spin stiffness, and the long-wavelength 
magnetization in the case of the susceptibility. The con- 
sistency among all the results obtained also contribute to 
our confidence in the procedures. 

Our final estimates for the critical couplings, taking 
statistically weighted averages of the values listed in Ta- 
bleU are g c = 2.52181(3) and 1.38882(2) for the full and 
incomplete bilayer, respectively. Thus we have improved 
the precision by two orders of magnitude relative to pre- 
vious calculations. Knowledge of the critical couplings to 
this level of accuracy should be useful for studies of vari- 
ous aspects of quantum criticality at low temperature in 
these systems, as one can avoid, to a higher degree than 
previously, effects from the eventual cross-over to renor- 
malixed classical or quantum-disordered behavior as de- 
viations from g c become relevant as T — > 0. 

For the correlation length, a weighted average of the 
four results from p s and Q2 listed in Table [I] gives 
v = 0.7106(9). This is consistent within error bars with 
the currently most accurate estimate of the 3D classi- 
cal Heisenberg exponent, v — 0.7112(5), obtained from 
classical 3D Heisenberg simulations in Ref. IT9I 

Taking a statistical average of the four individual re- 
sults for v (and analogously for g c discussed above) is mo- 
tivated here in spite of the fact that they are obtained in 
only two different simulations. This is because the wind- 
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ing numbers (giving p s ) and the sublattice magnetization 
(giving Q 2 ) are very weakly correlated in the simulations. 
We regard it as pure coincidence that the two values for 
each model in Tableware closer to each other than those 
of different models, and not an indication of a potentially 
different universality class in the incomplete bilayer (due 
to incomplete cancelation of Berry phasesiStlLSt) when 
the layer- exchange symmetry is not present. The very 
close agreement of the universal Binder ratio at g c also 
speaks in favor of the same universality class for both 
lattices. 

Although we have not quite reached the accuracy for 
v obtained in the most recent classical simulations^ (al- 
though our final error bar is actually only approximately 
twice as large), the precision is still sufficiently high to 
further increase the confidence in the belief that the uni- 
versality class of the transition is that of the 3D clas- 
sical Heisenberg model. The previously best (to our 



knowledge) determination of v for the transition in a 2D 
Heisenberg system is 0.70(1)^ 

We do not find the predicted^ universality of p s L 
at the critical point (whereas we do find the expected 
universality in the case of the Binder ratio); the values 
for the full and incomplete bilayer differ by about 30%. 
On the other hand, in recent simulations of the finitc- 
T transition of the 3D S — 1/2 Heisenberg ferromag- 
net and antiferromagnetj2& we find consistent values for 
both models at T c . Thus we are lead to speculate that 
there is some geometrical effect affecting the stiffness, so 
that universality might hold for different critical points 
(arrangement of couplings allowing a transition) on the 
same lattice but not necessarily for lattices with different 
unit cells. 
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